library(tidyverse)
library(patchwork)
library(emmeans)
library(simglm)
library(ggforce)
set.seed(3119)
sim_arguments <- list(
formula = y ~ 1 + hours + motivation + study + method,
fixed = list(hours = list(var_type = 'ordinal', levels = 0:15),
motivation = list(var_type = 'continuous', mean = 0, sd = 1),
study = list(var_type = 'factor',
levels = c('alone', 'others'),
prob = c(0.53, 0.47)),
method = list(var_type = 'factor',
levels = c('read', 'summarise', 'self-test'),
prob = c(0.3, 0.4, 0.3))),
error = list(variance = 20),
sample_size = 250,
reg_weights = c(0.6, 1.4, 1.5, 6, 6, 2)
)
df3 <- simulate_fixed(data = NULL, sim_arguments) %>%
simulate_error(sim_arguments) %>%
generate_response(sim_arguments)
test_study3 <- df3 %>%
dplyr::select(y, hours, motivation, study, method) %>%
mutate(
ID = paste("ID", 101:350, sep = ""),
score = round(y+abs(min(y))),
motivation = round(motivation, 2),
study = factor(study),
method = factor(method)
) %>%
dplyr::select(ID, score, hours, motivation, study, method)
study (two levels)test_study3 |>
group_by(study) |>
summarise(
avg_score = round(mean(score), 2)
)
mean_alone <- filter(test_study3, study == 'alone')$score |> mean()
mean_others <- filter(test_study3, study == 'others')$score |> mean()
study datatest_study3 |>
ggplot(aes(x = study, y = score, fill = study, colour = study)) +
geom_violin(alpha = 0.5) +
geom_sina(alpha = 0.5) +
theme(legend.position = 'none') +
stat_summary(colour = 'black', fun = mean, geom = 'point', size = 2) +
geom_segment(aes(x = 'alone', xend = 'others', y = mean_alone, yend = mean_others), colour = 'black') +
NULL
Plot on xy plane, dummy-coded.
contrasts(test_study3$study)
others
alone 0
others 1
alone = ref level.
xlim_lower <- -2.2
xlim_upper <- 2.2
ylim_lower <- -25
ylim_upper <- 55
test_study3 |>
mutate(study_num = ifelse(study == 'alone', 0, 1)) |>
ggplot(aes(x = study_num, y = score, fill = study, colour = study)) +
geom_jitter(alpha = 0.5, width = 0.1) +
scale_x_continuous(limits = c(xlim_lower, xlim_upper), expand = c(0, 0)) +
scale_y_continuous(limits = c(ylim_lower, ylim_upper), expand = c(0, 0)) +
geom_segment(aes(x = xlim_lower, xend = xlim_upper, y = 0, yend = 0), arrow = arrow(ends = 'both', length = unit(12, 'pt')), colour = 'black') +
geom_segment(aes(x = 0, xend = 0, y = ylim_lower, yend = ylim_upper), arrow = arrow(ends = 'both', length = unit(12, 'pt')), colour = 'black') +
geom_segment(aes(x = 0, xend = 1, y = mean_alone, yend = mean_others), colour = 'black') +
geom_segment(aes(x = 0, xend = 1, y = mean_alone, yend = mean_alone), colour = 'red', linewidth = 2) +
geom_segment(aes(x = 1, xend = 1, y = mean_alone, yend = mean_others), colour = 'red', linewidth = 2) +
NULL
Add note that the manual stuff is so you understand how coding works. But you’ll never have to do it this laboriously.
Exact same outcome:
contr.treatment()Aside: Foreshadow lecture 7, pm1 effect coding:
xlim_lower <- -2.2
xlim_upper <- 2.2
ylim_lower <- -25
ylim_upper <- 55
grand_mean <- mean(c(mean_alone, mean_others))
test_study3 |>
mutate(study_num = ifelse(study == 'alone', -1, 1)) |>
ggplot(aes(x = study_num, y = score, fill = study, colour = study)) +
geom_jitter(alpha = 0.5, width = 0.1) +
scale_x_continuous(limits = c(xlim_lower, xlim_upper), expand = c(0, 0)) +
scale_y_continuous(limits = c(ylim_lower, ylim_upper), expand = c(0, 0)) +
geom_segment(aes(x = xlim_lower, xend = xlim_upper, y = 0, yend = 0), arrow = arrow(ends = 'both', length = unit(12, 'pt')), colour = 'black') +
geom_segment(aes(x = 0, xend = 0, y = ylim_lower, yend = ylim_upper), arrow = arrow(ends = 'both', length = unit(12, 'pt')), colour = 'black') +
geom_segment(aes(x = -1, xend = 1, y = mean_alone, yend = mean_others), colour = 'black') +
geom_segment(aes(x = 0, xend = 1, y = grand_mean, yend = grand_mean), colour = 'red', linewidth = 2) +
geom_segment(aes(x = 1, xend = 1, y = grand_mean, yend = mean_others), colour = 'red', linewidth = 2) +
NULL
score ~ studymod1 <- lm(score ~ study, data = test_study3)
summary(mod1)
Call:
lm(formula = score ~ study, data = test_study3)
Residuals:
Min 1Q Median 3Q Max
-22.7902 -4.7902 0.2098 5.4766 18.2098
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 22.7902 0.6603 34.52 < 2e-16 ***
studyothers 4.7332 1.0093 4.69 4.53e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 7.896 on 248 degrees of freedom
Multiple R-squared: 0.08145, Adjusted R-squared: 0.07775
F-statistic: 21.99 on 1 and 248 DF, p-value: 4.525e-06
t.test(score ~ study, var.equal = T, data = test_study3)
Two Sample t-test
data: score by study
t = -4.6895, df = 248, p-value = 4.525e-06
alternative hypothesis: true difference in means between group alone and group others is not equal to 0
95 percent confidence interval:
-6.721058 -2.745251
sample estimates:
mean in group alone mean in group others
22.79021 27.52336
(mod1_emm <- emmeans(mod1, ~study))
study emmean SE df lower.CL upper.CL
alone 22.8 0.660 248 21.5 24.1
others 27.5 0.763 248 26.0 29.0
Confidence level used: 0.95
(p_mod1_emm <- mod1_emm |> plot() +
coord_flip())
mod1_emm_tib <- mod1_emm |>
as_tibble() |>
rename(score = emmean)
test_study3 |>
ggplot(aes(x = study, y = score, fill = study, colour = study)) +
geom_violin(alpha = 0.5) +
geom_jitter(alpha = 0.5, width = 0.2) +
theme(legend.position = 'none') +
geom_errorbar(data = mod1_emm_tib, aes(ymin = lower.CL, ymax = upper.CL), colour = 'black', width = 0.2) +
geom_point(data = mod1_emm_tib, colour = 'black', size = 3) +
NULL
levels(test_study3$study)
[1] "alone" "others"
For each comparison we want to make, assign 1 to the level we’re interested in, and assign -1 to the level you want to compare it to. If there are any other levels in there, assign them 0.
mod1_comparisons <- list(
'alone vs. others' = c(
1, # alone
-1 # others
)
)
contrast(
object = mod1_emm,
method = mod1_comparisons
)
contrast estimate SE df t.ratio p.value
alone vs. others -4.73 1.01 248 -4.690 <.0001
Basically reconstructs the model’s estimates. Not so interesting.
Original contrasts:
contrasts(test_study3$study)
others
alone 0
others 1
Now, let’s flip it:
contrasts(test_study3$study) <- contr.treatment(
levels(test_study3$study),
base = 2)
contrasts(test_study3$study)
alone
alone 1
others 0
xlim_lower <- -2.2
xlim_upper <- 2.2
ylim_lower <- -25
ylim_upper <- 55
test_study3 |>
mutate(study_num = ifelse(study == 'alone', 1, 0)) |>
ggplot(aes(x = study_num, y = score, fill = study, colour = study)) +
geom_jitter(alpha = 0.5, width = 0.1) +
scale_x_continuous(limits = c(xlim_lower, xlim_upper), expand = c(0, 0)) +
scale_y_continuous(limits = c(ylim_lower, ylim_upper), expand = c(0, 0)) +
geom_segment(aes(x = xlim_lower, xend = xlim_upper, y = 0, yend = 0), arrow = arrow(ends = 'both', length = unit(12, 'pt')), colour = 'black') +
geom_segment(aes(x = 0, xend = 0, y = ylim_lower, yend = ylim_upper), arrow = arrow(ends = 'both', length = unit(12, 'pt')), colour = 'black') +
geom_segment(aes(x = 0, xend = 1, y = mean_others, yend = mean_alone), colour = 'black') +
geom_segment(aes(x = 0, xend = 1, y = mean_others, yend = mean_others), colour = 'red', linewidth = 2) +
geom_segment(aes(x = 1, xend = 1, y = mean_others, yend = mean_alone), colour = 'red', linewidth = 2) +
NULL
mod1b <- lm(score ~ study, data = test_study3)
coef(mod1b)
(Intercept) studyalone
27.523364 -4.733155
Compare this to the original coefs from mod1:
coef(mod1)
(Intercept) studyothers
22.790210 4.733155
(mod1b_emm <- emmeans(mod1b, ~study))
study emmean SE df lower.CL upper.CL
alone 22.8 0.660 248 21.5 24.1
others 27.5 0.763 248 26.0 29.0
Confidence level used: 0.95
p_mod1b_emm <- mod1b_emm |> plot() +
coord_flip() +
ggtitle('With reference level = others')
p_mod1_emm <- p_mod1_emm + ggtitle('With reference level = alone')
p_mod1b_emm + p_mod1_emm
Identical!
method (three levels)test_study3 |>
group_by(method) |>
summarise(
avg_score = round(mean(score), 2)
)
mean_read <- filter(test_study3, method == 'read')$score |> mean()
mean_self <- filter(test_study3, method == 'self-test')$score |> mean()
mean_summ <- filter(test_study3, method == 'summarise')$score |> mean()
test_study3 |>
ggplot(aes(x = method, y = score, fill = method, colour = method)) +
geom_violin(alpha = 0.5) +
geom_jitter(alpha = 0.5, width = 0.2) +
theme(legend.position = 'none') +
NULL
Goal: Compare the means of each group to one another. That was easy when we had only two groups: one straight line can go from mean to mean. But now one straight line won’t cut it anymore.
The smallest number of lines that will connect all the means = 2. Instead of estimating one line, we will estimate two.
And in general, if we have \(n\) groups, we will be estimating \(n-1\) lines.
The lines that we’re going to be estimating now:
test_study3 |>
ggplot(aes(x = method, y = score, fill = method, colour = method)) +
geom_violin(alpha = 0.5) +
geom_jitter(alpha = 0.5, width = 0.2) +
theme(legend.position = 'none') +
stat_summary(colour = 'black', fun = mean, geom = 'point', size = 2) +
# line from read to self-test:
geom_segment(colour = 'black', aes(x = 'read', xend = 'self-test', y = mean_read, yend = mean_self)) +
# line from read to summarise:
geom_segment(colour = 'black', aes(x = 'read', xend = 'summarise', y = mean_read, yend = mean_summ)) +
NULL
But I’m lying a little bit. This looks like going to go from 0 to 1 in one line, and 0 to 2 in the other line. But that’s not true – this is just the intuition about what groups we are comparing. Really what we’re going to do is this:
method in xy spaceTwo panels: the first will be the difference between
read and self-test (the first variable), the
second will be difference between read and
summarise (the second variable).
xlim_lower <- -2.2
xlim_upper <- 2.2
ylim_lower <- -25
ylim_upper <- 55
p1 <- test_study3 |>
filter(method %in% c('read', 'self-test')) |>
mutate(method_num = ifelse(method == 'read', 0, 1)) |>
ggplot(aes(x = method_num, y = score, fill = method, colour = method)) +
geom_jitter(alpha = 0.5, width = 0.1) +
scale_x_continuous(limits = c(xlim_lower, xlim_upper), expand = c(0, 0)) +
scale_y_continuous(limits = c(ylim_lower, ylim_upper), expand = c(0, 0)) +
geom_segment(aes(x = xlim_lower, xend = xlim_upper, y = 0, yend = 0), arrow = arrow(ends = 'both', length = unit(12, 'pt')), colour = 'black') +
geom_segment(aes(x = 0, xend = 0, y = ylim_lower, yend = ylim_upper), arrow = arrow(ends = 'both', length = unit(12, 'pt')), colour = 'black') +
geom_segment(aes(x = 0, xend = 1, y = mean_read, yend = mean_self), colour = 'black') +
geom_segment(aes(x = 0, xend = 1, y = mean_read, yend = mean_read), colour = 'red', linewidth = 2) +
geom_segment(aes(x = 1, xend = 1, y = mean_read, yend = mean_self), colour = 'red', linewidth = 2) +
labs(
title = 'First line: read vs. self-test'
) +
NULL
p2 <- test_study3 |>
filter(method %in% c('read', 'summarise')) |>
mutate(method_num = ifelse(method == 'read', 0, 1)) |>
ggplot(aes(x = method_num, y = score, fill = method, colour = method)) +
geom_jitter(alpha = 0.5, width = 0.1) +
scale_x_continuous(limits = c(xlim_lower, xlim_upper), expand = c(0, 0)) +
scale_y_continuous(limits = c(ylim_lower, ylim_upper), expand = c(0, 0)) +
geom_segment(aes(x = xlim_lower, xend = xlim_upper, y = 0, yend = 0), arrow = arrow(ends = 'both', length = unit(12, 'pt')), colour = 'black') +
geom_segment(aes(x = 0, xend = 0, y = ylim_lower, yend = ylim_upper), arrow = arrow(ends = 'both', length = unit(12, 'pt')), colour = 'black') +
geom_segment(aes(x = 0, xend = 1, y = mean_read, yend = mean_summ), colour = 'black') +
geom_segment(aes(x = 0, xend = 1, y = mean_read, yend = mean_read), colour = 'red', linewidth = 2) +
geom_segment(aes(x = 1, xend = 1, y = mean_read, yend = mean_summ), colour = 'red', linewidth = 2) +
labs(
title = 'Second line: read vs. summarise'
) +
NULL
p1 + p2
method variableManual dummy-coding:
dummyCoded <- test_study3 %>%
select(ID, score, method) %>%
mutate(
method1 = ifelse(method == "self-test", 1, 0),
method2 = ifelse(method == "summarise", 1, 0))
dummyCoded
contrasts(test_study3$method)
self-test summarise
read 0 0
self-test 1 0
summarise 0 1
score ~ methodmod2 <- lm(score ~ method, data = dummyCoded)
summary(mod2)
Call:
lm(formula = score ~ method, data = dummyCoded)
Residuals:
Min 1Q Median 3Q Max
-23.4138 -5.3593 -0.1959 5.7496 17.8041
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 23.4138 0.8662 27.031 <2e-16 ***
methodself-test 4.1620 1.3188 3.156 0.0018 **
methodsummarise 0.7821 1.1930 0.656 0.5127
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 8.079 on 247 degrees of freedom
Multiple R-squared: 0.04224, Adjusted R-squared: 0.03448
F-statistic: 5.447 on 2 and 247 DF, p-value: 0.004845
(Nearly) equivalent t-tests (different bc diff degrees of freedom, bc filtering dataset. but close).
dC1 <- test_study3 |>
filter(method %in% c('read', 'self-test')) |>
mutate(method = factor(method, levels = c('read', 'self-test')))
contrasts(dC1$method) <- contr.treatment(2)
contrasts(dC1$method)
2
read 0
self-test 1
t.test(score ~ method, var.equal = T, data = dC1)
Two Sample t-test
data: score by method
t = -3.1354, df = 151, p-value = 0.002063
alternative hypothesis: true difference in means between group read and group self-test is not equal to 0
95 percent confidence interval:
-6.784657 -1.539272
sample estimates:
mean in group read mean in group self-test
23.41379 27.57576
dC2 <- test_study3 |>
filter(method %in% c('read', 'summarise')) |>
mutate(method = factor(method, levels = c('read', 'summarise')))
contrasts(dC2$method) <- contr.treatment(2)
contrasts(dC2$method)
2
read 0
summarise 1
t.test(score ~ method, var.equal = T, data = dC2)
Two Sample t-test
data: score by method
t = -0.66342, df = 182, p-value = 0.5079
alternative hypothesis: true difference in means between group read and group summarise is not equal to 0
95 percent confidence interval:
-3.108082 1.543915
sample estimates:
mean in group read mean in group summarise
23.41379 24.19588
Before, I said that the smallest number of lines we need to compare three groups is two. But this actually means that one of the possible comparisons is left out.
(??? above the missing line between self-test and summarise)
Options:
(mod2_emm <- emmeans(mod2, ~method))
method emmean SE df lower.CL upper.CL
read 23.4 0.866 247 21.7 25.1
self-test 27.6 0.994 247 25.6 29.5
summarise 24.2 0.820 247 22.6 25.8
Confidence level used: 0.95
mod2_emm |> plot() +
coord_flip()
mod2_emm_tib <- mod2_emm |>
as_tibble() |>
rename(score = emmean)
test_study3 |>
ggplot(aes(x = method, y = score, fill = method, colour = method)) +
geom_violin(alpha = 0.5) +
geom_jitter(alpha = 0.5, width = 0.2) +
theme(legend.position = 'none') +
geom_errorbar(data = mod2_emm_tib, aes(ymin = lower.CL, ymax = upper.CL), colour = 'black', width = 0.2) +
geom_point(data = mod2_emm_tib, colour = 'black', size = 3) +
NULL
levels(test_study3$method)
[1] "read" "self-test" "summarise"
Assign 1 to the level we are interested in, and -1 to the baseline we want to compare it to.
mod2_comparisons <- list(
'Self-test vs. Read' = c(-1, 1, 0), # read self-test summarise, in that order
'Summarise vs. Read' = c(-1, 0, 1),
'Self-test vs. Summarise' = c(0, 1, -1),
'mean(Self-test, summarise) vs. Read' = c(-1, 0.5, 0.5) # weights must sum to 0 -- pools self-test and summarise tog and pits that mean against mean of Read
)
1 - 1
mean(Self-test) - mean(Read) = positive, bc self-test is
larger than read
(mod2_comparisons_test <- contrast(mod2_emm, mod2_comparisons))
contrast estimate SE df t.ratio p.value
Self-test vs. Read 4.162 1.32 247 3.156 0.0018
Summarise vs. Read 0.782 1.19 247 0.656 0.5127
Self-test vs. Summarise 3.380 1.29 247 2.622 0.0093
mean(Self-test, summarise) vs. Read 2.472 1.08 247 2.290 0.0229
contrast() gives us the p-values for each difference. We
can also get the 95% CIs of each difference by giving the outcome of
contrast() into confint().
confint(mod2_comparisons_test)
contrast estimate SE df lower.CL upper.CL
Self-test vs. Read 4.162 1.32 247 1.564 6.76
Summarise vs. Read 0.782 1.19 247 -1.568 3.13
Self-test vs. Summarise 3.380 1.29 247 0.841 5.92
mean(Self-test, summarise) vs. Read 2.472 1.08 247 0.345 4.60
Confidence level used: 0.95